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Point-ion models have been extensively used to determine "hole numbers" at copper and oxygen 
sites in high-temperature superconducting cuprate compounds from measured nuclear quadrupole 
frequencies. The present study assesses the reliability of point-ion models to predict electric field 
gradients accurately and also the implicit assumption that the values can be calculated from the 
"holes" and not the total electronic structure. First-principles cluster calculations using basis sets 
centred on the nuclei have enabled the determination of the charge and spin density distribution 
in the Cu02-plane. The contributions to the electric field gradients and the magnetic hyperfine 
couplings are analysed in detail. In particular they are partitioned into regions in an attempt to 
find a correlation with the most commonly used point-ion model, the Sternheimer equation which 
depends on the two parameters R and 7. Our most optimistic objective was to find expressions for 
these parameters, which would improve our understanding of them, but although estimates of the 
R parameter were encouraging the method used to obtain the 7 parameter indicate that the two 
parameters may not be independent. The problem seems to stem from the covalently bonded nature 
of the Cu02-planes in these structures which severely questions using the Sternheimer equation for 
such crystals, since its derivation is heavily reliant on the application of perturbation theory to 
predominantly ionic structures. Furthermore it is shown that the complementary contributions of 
electrons and holes in an isolated ion cannot be applied to estimates of electric field gradients at 
copper and oxygen nuclei in cuprates. 

I. INTRODUCTION 

There is a large quantity of nuclear magnetic resonance (NMR) data from high-temperature superconducting 
cuprate crystals from which electric field gradients (EFG) can be derived. EFG's are a measure of the non-spherical 
components of the charge distribution surroundi™g-|tJ|ic of interest and is used to estimate the hole population 

in models of superconductivity. Most estimationsoHijIj'QJuyul have been made using a point-ion model with Sternheimer 
correction factors, called here the Sternheimer equationBcl {SE^l which is. briefly discussed in Sec. II. In particular the 
measured changes of the EFG on doping have been discussedli3EilOllj'Eil in terms of the distribution of the additional 
holes among_Lhe orbitals on each ion. In related areas of research, however, point-charge models are apparently no 
longer in uselij although there does not seem to be any report in the literature discussing the unreliability of the SE. In 
this paper we address this problem, particularly for systems where there is evidence that the bonding between "ions" 
is more covalent than ionic. The Cu02 sheets, which are a common feature of all cuprates showing high-temperature 
superconductivity behaviour, are thought to have bonds showing a distinct covalent character. 

In addition to the point-ion approximation, the above mentioned semi-empirical analyses of EFG values in cuprate 
superconductors are also based on the assumption that the EFG values at a nucleus can be calculated (with opposite 
sign) for that configuration where the unoccupied spin orbitals are assumed to be occupied and the occupied spin 
orbitals are assumed to be unoccupied. For the Cu ion, in particular, it is expected that the EFG is just the 
same as for but with opposite sign. It seems that this concept is widely adopted unconditionally. In the present 
paper we also address this assumption, called here the electron-hole symmetry. We demonstrate that it is entirely 
unjustified for the Cu ions in a cuprate environment but also leads to false estimates of the EFG values at the oxygens 
in the Cu02-planes. 

In order to quantify our doubts on the applicability of the SE and of the electron-hole symmetry to the evaluation 
of EFG in copper oxides we have used the wave functions from previously published first-principles calculat ions 
and complemented them for illustrative purposes with additional simulations. These are cluster calculationsc£ltZI 
on La2Cu04 and YBa2Cu307 cuprates using the density functional theory with local density approximation and 
generalized gradient corrections, which have provided EFG data for the Cu and O in the GiLQa-planes in agreement 
with experiment. Calculations which used augmented plane waves give similar agreementll§E3 for all nuclei except 
copper in the Cu02 plane. It is not the purpose of the present paper to discuss such differences since the focus is 
solely an-jthe reliability of point-ion calculations. For detailed comparisons with experimental measurements we refer 
to Refsm. 

The general idea of the cluster approach to electronic structure calculations of properties which depend upon 
predominantly local electron densities is that the parameters that characterize a small cluster should be transferable 
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to the solid and largely determine its properties. The essential contributions to EFGs and to magnetic hyperfine fields 
are given by rather localized interactions and therefore it is expected that these local properties can be determined 
and understood with clusters calculations. Approximations must be made concerning the treatment of the lattice 
in which the cluster is embedded. Using as large a cluster as is possible is of course advantageous. It is necessary, 
however, that the results obtained should be checked with respect to their depe ndence on the cluster size. 

and the general contributions to the 



The basic principles of cluster calculations are briefly discussed in Sec. [II A 



one-electron operator from regional partitioning are given in Sec. [II B 
operator and also to the hyperfine coupling operator 



In Sec. 



[II C, this is then apphed to the EFG 



The latter which can be used for clusters with unpaired electron 
spins is very similar to the EFG operator. The only difference is that it uses the spin density and the EFG the charge 
density. 

Effectively the cluster calculation of the EFG is divided into contributions from the ion of interest (the target 
ion), the rest of the cluster and the overlap between these two. This regional partitioning technique is described in 
Sec. 1^ and the contributions to the EFG for a particular cluster calculation are given as an example. In Sec. ^ 
correlations of first-principles partitions with Sternheimer terms are investigated. The above mentioned partition 
enable th e Ste rnheimer antishielding factors, R and 7 to be associated with quantities calculated from first principles 
(see Sec. |V B| ), which in turn allows us to compare the predictions of the SE for a model cluster with that of the ab 
initio calculation of t he same cluster. This provides a much more sensitive test of the SE than could otherwise be 
obtained. In Sec. V C the contributions from the target ion are analysed in terms of the individual orbitals indicating 
how the "holes" 

In Sec. 



VD and VE 



have been determined. The shortcomings of the simplifying approaches are pointed out. It is shown 
that the values of the EFG are determined by a subtle cancellation of large individual terms, 
the electron-hole symmetry is studied. A summary and conclusions are given in Sec. VI. 

Except for energies, atomic units are used throughout, i.e. the EFG components Vu are given in e Ug'^ 
— Ha/a^. The quantity qu ~ Viij\e\ then corresponds to —9.7174 x 10^^ in imits of V/m^. 
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II. THE STERNHEIMER EQUATION 



The Sternheimer equation has been written in the following for: 




y^^ = (1 - R)vl°'^''^ + (1 - (2.1) 

where Ya is one of the diagonal components of the EFG tensor for a target ion which can be determined experimentally, 
y^ocai experimental EFG component of the target free ion and ]/A'^***'=«= is the contribution to the EFG component 

from the charges in the lattice surrounding the target ion. The two parameters in the equation are both antishielding 
factors; R arises from the electrons in the valence shell of the target ion which arc possibly overlapping the electron 
distribution of the nearest neighbour ions, and 7 accounts for the contribution of the EFG due to the polarisation of 
the target ion in the electric field of the environmental charges. 

The crystal structure is therefore split up into three regions, the isolated target ion (here referred to as local), the 
space between the target ion and its nearest neighbours, the rest of the ions in the crystal (here referred to as lattice). 

There is large spread of values for the parameters, in particular 7, derived from EFG df||fhft'r superconducting 
cuprates. For copper, values for 7 of -7.6, -10.4, -17, and -20 have been reported in Refs.liH3'BQ, respectivelsz-. For 
planar oxygen ions the lattice contributions (1 — 7)Vj"""^'^ accounts for 36% and 60% of the total Vxx hi Ref.E3 and 
Ref.EJ, respectively. 



III. FIRST PRINCIPLES CALCULATIONS ON CLUSTERS 



A. Description of cluster calculations 



A cluster is a careful selection of ions within a crystal which are intended to be able to calculate localised properties 
accurately. The target ion (the ion whose properties are to be calculated) should be at or very near the centre of 
the cluster. The target and normally, at least, its nearest neighbours form the core of the cluster and are treated 
most accurately using first-principles all-electron methods. Outside this core the next shell of positively charged 
ions are represented by pseudopotential functions which have been shown to behave better than just bare charges to 
represent the ions since pseudopotentials prevent unrealistic electron density distortions characteristic of the positive 
point charges. Bare charges (~ 2000) are used outside the shell of pseudopotentials to simulate the rest of the crystal 
lattice. Some of the more remote charges from the target ion are moved slightly so that the target ion experiences 
the correct Madelung potential. 



3 



The present work aims at an assessment of calculations of EFG with a point-ion model and Sternheimer correc- 
tions and comparison with first-principles methods. To illustrate the problems we use here the results from three 
different clusters which all simulate the compound La2Cu04 but we note that similar results have been obtained for 
YBa2Cu307. The core of these clusters comprise one, two and nine copper ions, respectively, each with an appropriate 
number of nearest neighbour oxygen ions. 

Only the central Cu ion in the cluster Cu06/Cu4Laio (Fig. |l|(a) where X=Y=La) is used as the target ion which, 
together with the 6 nearest neighbour oxygen ions, forms the core of the cluster used for the all-electron calculation. 
The neighbouring 4 Cu and 10 La ions are represented by pseudo-potential functions. In the Cu20ii/Cu6Lai6 cluster 
(Fig. 0(b) where X=Y==Z=La) we have used both the central oxygen ion and the neighbouring two Cu ions separately 
as target ions. The core of the cluster additionally includes the 10 nearest neighbour oxygen ions. The adjacent 6 Cu 
and 16 La ions are represented by pseudo-potential functions. 

The nuclear positions have been chosenE3 according to the tetragonal structure of La2Cu04 (space group lA/mmm) 
with a = b = 3.77 A and c — 13.18 A, with a Cu-O(a) distance of 2.40 A and with a Cu-La distance of 4.77 A. 

The cluster core uses a 6-311G basis set as provided by Gaussian 98. The density functional method was used to 
obtain the wave functions from which the EFG at the target ions were calculated. This procedure has been used 
consistently by us since the wave functions which are produced give calculated properties which are in agreement with 
experimental values. 

It should be noted that the lattice region of the Sternheimer equation includes all bare charges, the pseudo-potential 
ions and all the ions of the cluster core except the target ion. 



B. Contributions to the one-electron operator from regional partitioning 

Let us consider a system of N nuclear centres. The X*'' centre, at site Rk, is the origin for tik basis functions 
which are mutually orthogonal. The /c*'' basis function on site Rk is denoted by BK,k{r — Rk)- The total number of 
basis functions (atomic orbitals) is 

TV 

J^c = ^ nK- (3.1) 

K=l 

the c in ric identifies that our system here is the core of a cluster. The molecular orbitals (MO), (/>, of the system 
are orthogonal linear combinations of the atomic orbitals. We allow for two sets of MO's, one set to hold electrons of 
a-spin projection, the other set /3-spin projection. The m}^ MO of a-spin projection is 

N N 7lK 

cl^mA^ = E Rk) = J2J2''^i^KAr~ Rk) (3.2) 

K=l K=lk=l 

where the c's are the MO coefficients. 

The expectation value of any quantity, corresponding to the operator 0{f), associated with a nuclear site Rj for 
the MO 4>m,a{'r)i is given by the matrix element 

MZ,o.{Rj) =< 4>r,^A^\0{r- Rj)\4>,nA^) > ■ (3.3) 
Developing the MO's according to Eq. ( ^ ) we get 

M^ARj) -EE r™,a(^, L) (3.4) 

K L 

where, for convenience, we have defined: 

tZAK,L) = EE^m^^'a < BkAt- RK)\0{r- Rj)\Bl,{t- Rl) > ■ (3.5) 

k I 



We note that in Eq. (3.4) K and L sum over all nuclear centres of the cluster and that the target ion is J. The SE 
concentrates entirely on the target ion and so will we. Hence in Eq. (^^) we separate out the target ion as follows: 

M^JRj) = ^ZAJ,J)+J2 r™,a(^>>/) 

+ Er™,a(^,i) + E Er™,-(^'^) (3.6) 
= 'mZARj) + '^mZARj) + '''mZJRj) (3-7) 
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noting that ''M^jRj) = Ek^.j ^ZJK, J) + Y^l^j r«,„( J, L). 
This results in the identification of three diff'erent terms. 

1. The first term comprises ah contributions from on-site basis functions (that is, ah basis functions centred at the 
target ion Rj) and is denoted by regional partition I. 



2. The second and third terms in Eq. ( |3.6D are numerically identical and contain contributions arising from both 
on-si te a nd off-site [Rk-K ^ J) basis functions (corresponding to regional partition II) and denoted by II in 
Eq. (13). 

3. ^^^Af^ Q,(_Rj) contains no reference to the on-site basis functions. 

This is illustrated schematically in Fig. ^. 
If O is the identity operator 1 then 

E E L) = E E E ^^'X^nt < BkA^- RK)\BLAr- Rl) > 

K L K,L k I 

= 'MijRj) + ''MijRj) + '''M,ijRj) (3.8) 

and since the basis orbitals (functions) on each centre have been conveniently chosen to be orthogonal 

"J 

'Mi^Rj) = E(c;!;'a)' = 'N^,M)- (3-9) 

k 

If we define an overlap integral as SK,k,L,i BK,k{r — RK)\BL.i{r — Rl) > we get 

''MijRj) = 2 E E ^-^'c'^^iSj^k^Kj = '%nAJ)- (3-10) 

K^J k 

Whereas ^Nm,a{J) has been interpreted by Mulliken as the charge on atom J due to a-spin electrons in MO (/)m,a, 
^^Njn a{J) is the a-spin density that atom J shares with all its neighbours in the same MO. With the definition 



Pa{Rj) E Pm.ARj) = E 



occ r 



(3.11) 



(noting that the sum is only over occupied MOs only) this gives the Mulliken charge density attributed to the nuclear 
centre at Rj: 

pMull^jl^^ ^ p^(Rj) + p^{Rj), (3.12) 

It should be noted that the Mulliken analysis of the charge distribution has no physical meaning but it is very 
useful when discussing the charge distribution in molecules and clusters, perhaps more suited to systems which are 
non-ionic rather than ionic. 



C. The EFG and hyperfine coupling operators 

The concepts developed above are now applied to expectation values of the operator 



V'^x) = (V,Vj - h,A)- - %A-. (3.13) 



The expectation values of this operator cover three contributions of interest: 
1. The Fermi contact density. 



This is just the second term on the right hand side of Eq. ( 3.13 ) and its expectation value gives rise to the 



expression ^1-0(0)1^, where 1-0(0)1^ is the spin density at the target nucleus. We will not discuss this term in 
the following but we note that the same analysis which will be performed for the contributions to the EFG has 
also been applied to the contact density. The corresponding results are given in Appendix |^ 
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2. EFG operator 



This operator is only the first term in Eq. ( 3.13 ) and transforms as a spherical harmonic of order 2 and its 



expectation value for an s-like charge or spin distribution vanishes. It is written as 

P^^(.f) = ^^^^i-^^. (3.14) 

The expectation values of 

V'i =V'^{x- R.j) (3.15) 

determine the EFG tensor only if the total charge density distribution is used. 

3. Dipolar Hyperfine Coupling operator 

This has exactly the same form as the EFG operator, but is only used with the total spin density distribution. 

It will be convenient, and unlikely to cause confusion since the only other operator defined is the unit operator J, 



to replace O in equations ( p.3| - 3.7) with just ij when we should write V j (there is no need to repeat J) 



= h'P„\ARj) + "M^, Ji?j) + "^M^,,„(i?j). (3.16) 

Since the operator (15) contains a factor roughly proportional to the reciprocal of the cube of the distance from 
nuclear centre J, we would expect only those terms which describe the electron density close to the centre J would 
be significant. Clearly one of these terms is ^M^l ^.(-Rj). Explicitly 

n 7 n 7 

'M%^M -Y.T.'''^,yd^DlMiiJ) (3-17) 

k I 

where D^k k l li'^) Bx.ki^ ^ Rk)\'^'^j \Bl.i{x — Rl) >■ By summing over all the occupied MO's we define the 
quantities 

occ occ 

^Gt = E ^^it^^^R-^) 'G^p = E 'M^^i^j)- (3-18) 

m m 

The contribution of the "on-site" terms (I) to the EFG is then given by the sum 
and the difference 

is the corresponding contribution to the dipolar hyperfine tensor T. 

Analogous definitions determine the mixed on-site and off-site contributions (II) and those of purely off-site contri- 
butions (III). This leads to the following representations for the total EFG tensor 

Y^j ^ ly^] llY^j niyu _^ T^^i (■3 21) 

and for the dipolar term 

rjiij _ Irpij ^ Urpij ^ Hl/jiij (3 22) 



Note that the last term in Eq. ( 3.21 ) represents the EFG contribution of all nuclear point charges Zk- It is given 

by 

W'^ = ^ ^C^-g ~ Rj)i{RK - Rj)] - 5ij\RK - -Rjp (3 23) 

K^j \Rk - RjY' 

Of course W^^ has no place in the electron-nuclear hyperfine coupling tensor. 
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D. Density Matrix Formulation of Partitioning 



The previous two subsections can be more succinctly described using density matrix terminology. In fact in Sec. VE 
it enables certain conclusions to be reached which would be difficult to achieve otherwise. 

Let B be the column matrix of the complete set of basis functions Bfc i^{f— Rk)- It will be made clear shortly why 
we want to order the basis functions such that those belonging to the target ion are placed at the beginning. The 
MO's of a-spin projection can be written as the column matrix 



(3.24) 



where Cq. is the matrix of the MO coefficients, each column corresponding to a particular MO. Equation (3.2) corre- 
sponds to the m*'* row of $q- which can be written here as 4>m,a = Cm where Cm,a is the m*'* column of Cq. 
In order to obtain the expectation values of the operator 0{f) it is convenient to define the matrix 



b*^ B|0(r)|Bt > 



(3.25) 



so that 



(3.26) 



which is equivalent to Eq. ( |3.3| ). Pm is the density matrix Crn.ac\n a- accordance with previous practice we have 
similarly for electrons with /3-spin projection Q,„ = Cm^^cJ^ ^. 

Instead of discussing just individual MO's we can usefully define a total density matrix for all the electrons with 
a-spin projection as 



r T 



(3.27) 



where Iq is a diagonal matrix with 1 on the m*'* diagonal if MO m is occupied and zero otherwise. 

Since we have ordered the ric basis orbitals of the cluster core such that all the rit basis orbitals on the target ion 
are listed first we partition the density matrix as follows, 



III 



(3.28) 



Pm / 



where ^Pm is an rit x rit matrix, ^^Pm. is an fit x {ric — rit) matrix and ^^^Pm is an (ric — rit) x {ric — nt) matrix. If we 
define 



\ lie — nt ^nt 



ntOnc-nt 



c — ntOric-nt / 



(3.29) 



Hp _ 



\ Pm ■ iic-ntOnc-nt / 



(3.30) 



and 



nil 



( 

\ nc-ntOnt 



nt^ ric-nt 



III 



Prn / 



(3.31) 
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where jOk is a null j x k matrix, we have 

Pm ^ Pru ^ Pm ^ Pm- (3.32) 

Hence we have an equivalent expression for Eq. ( ^.7| ) 

M^ciRj) = Tr(ip„b'^) + Tr("P„b'^) + rr("iP™b°) (3.33) 

that is, ^M^ jRj) = Tr(Rp„b'^). 

Having demonstrated the equivalence between the two mathematical approaches, for example, Eq. ( 3.1^ ) can be 
rewritten equivalently as 

^G'^ ^Tri^PW^) ^G'^ ^Tri^QW^). (3.34) 

Thus the density matrix approach emphasises the partition method chosen diagrammatically. 
Eqns. (3. IE) and ( 3.20] ) can be succinctly written as 



=rr((ip+iQ)b*^) (3.35) 

and 

lyu = Tr((ip -iQ)b*^) (3.36) 
where (^P +^Q) is the total charge density matrix and (^P — ^Q) is the total spin density matrix for region I. 

IV. AN EXAMPLE OF REGIONAL PARTITIONING 

To investigate the electronic structure of La2Cu04v a parent compound of high-temperature superconducting ma- 
terials such as La2-xSra;Cu04, we have performedliB extended first-principles cluster calculations. Several clusters 
containing up to nine copper atoms embedded in a background potential were investigated. In Fig. ^ the highest 
occupied molecular orbital of the cluster Cug042/Cui2La5o is shown. All electron triple-zeta basis sets (6-311G 
Gaussian functions) were used for nine Cu and 42 O atoms resulting in a total of 663 electrons. 

The detailed results of the spin polarized calculatio ns us ing the local cluster approximation with generalized gradient 
corrections (BLYP functional) will be given in Sec. VC. Anticipating these results, here the contributions to 



and from the three regional partitions R = I, H, and HI, respectively, for the central copper and the oxygen 

atom indicated in Fig. || are collected in Table |. 

In contrast to point-ion charge models where only the (small) valence charge is considered, the electronic structure 
is here determined by using all-electron basis sets (including core electrons) on the atoms in the centre of the cluster. 
The contributions of the nuclear point charges W is, however, cancelled to a large extent by the off-site contributions 
^^^Ga.fs- Adding and subtracting the contributions from the a and P spin projections we get the EFG components and 
dipolar hyperfine couplings, respectively, listed in Table 0. It is seen that the combined contributions from region III 
and W are small. For the oxygen, the values from region II give a reduction of the main contributions from region I 
by about 10 %, for copper by 20 %. For the dipolar hyperfine couplings, the contributions from regions other than I 
mostly cancel. 

It is remarkable that the calculations on the small clusters shown in Fig. |l| already give values for the EFG and the 
dipolar hyperfine couplings that are close to those obtained from the large cluster with nine copper ions. With the 
cluster Cu06/Cu4Laio (Fig. |^(a) where X=Y=La) we obtain at the copper the values Vzz = 1.396 and T^^ — —3.526. 
With the cluster Cu20ii/Cu6Lai6 (Fig. |l|(b) where X=Y=Z=La) we get V^z = 1.167 and T^z = -3.467 and, for the 
O as target ion, Vxx = —0.862 and T^x — 0.652. This demonstrates that these properties depend on the local charge 
and spin distributions and that cluster approaches are especially suited for their detailed investigations. 



V. CORRELATION OF FIRST-PRINCIPLES PARTITIONS WITH STERNHEIMER TERMS 



A. Introductory remarks 



Ideally we would like to have a correspondence between the terms of the first-principles calculation, Eq. (3.21), 
with the semi-empirical Sternheimer equation, Eq. (2.1). But the approaches are quite different. The first-principles 
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approach is a straightforward apphcation of molecular orbital theory using the density fun ction al method to obtain 
the electron density from which the expectation value of the appropriate operator (Eq. ( |3.13| )) is calculated. We 
would like to emphasise that the particular choice of theoretical method is not crucial. The same argument applies 
for Hartree-Fock or improved methods like multi-configuration self-consistent field, M0ller-Plesset and configuration 
interaction methods. The limitations of the approach are largely determined by available computer resources which 
in our case effectively determine the size of the cluster we can use. 

In the Sternheimer approach the starting point is to regard the target ion in a crystal as isolated and then add terms 
to compensate for the recognised interactions. Although historically the semi-empirical approach always precedes the 
first-principles methods, sometimes by many years, the value of the Sternheimer equation should not be discarded 
lightly since it appeals to the perturbed atomic picture, a model which has been and still is the cornerstone of 
experimental chemistry, rather than the more intractable molecular picture. On the other hand, as in chemistry 
where some molecules demonstrate high degrees of electron delocalisation, it is necessary to abandon the atomic 
picture in favour of special molecular models. A good example is the separation of aliphatic and aromatic organic 
chemistries. However the subject matter is enormous in both these areas and the separation is more than justified. 
This is to be contrasted with the relatively small number of high-Tc superconducting materials so it was natural to 
pursue the perturbed atom approach. It is clear that because of the above differences in approach an immediate 
correspondence between the first principles approach and the semi-empirical approach is not to be found. But since 
we understand the terms in the SE there is a chance of some correspondence. 



B. Identification of terms in the Sternheimer Equation 



Firstly we look a t the VA''***'=<^ term in Eq. (2.1). This is the contribution to the EFG of the environment of the 
target ion. In Eq. (3.21) this is the purely off-site contribution of the cluster, ™Vii, and the contribution from the 



point charges surrounding the cluster, Wa (Eq. ( |3.23| )): 

ylatUce ^Illy^^^^^^^ (5.1) 

yiocal atomic (ionic) term which is at the core of the semi-empirical perturbation approach. Cluster calculations 
are analogous to molecular orbital calculations where the properties of atoms largely disappear as identifiable entities 
although some analyses of th£_jelectron density distribution attempt to allocate charges to particular atoms (such as 
the commonly used MuUikencj population analysis). However it must be stressed that such analyses serve only as 
a useful guide, and are not without controwjrsy. One studyE3 concLtided that the LowdinEj population analysis was 
more appropriate than either the MuUikcnEil or a modified Mullikenc3 designed for heteronuclear bonds. 

Previous cluster calculations have shown that the degree of overlap and delocalisation of the target ion electrons is 
considerable which places in doubt the validity of perturbation methods using an isolated ion as the zeroth function. 

Vii is the closest we can get to since Vu is the contribution to Va from only the basis orbitals centred on 

the target atom. 

The semi-empirical approach is well aware of the potential distortion, and associated anti-shielding effect on Va, 
of the target ion by the environment of charges. The R and 7 terms were introduced to accommodate this. Since 
the first-principles calculation purports to include such distortions automatically in ^Vu we are inclined to absorb 
—jyio-tttce j^^Q yAoca' replace it in Eq. (pTl) by ^Vu to give approximately: 



Vu = {1-R) 'Vu + Wu. (5.2) 



However comparing Eq. ( |5.2| ) with Eq. (3.21), noting that we have already identified the first-principles correlation 
of Vl^'^*""', we deduce that 

-R'V, = 'V,, (5.3) 

enabling us to estimate R from our first-principles regional partitioning approach. We obtain (values used in some 
semi-empirical estimates are in brackets) 

Rcu = 0.21 (0.2) and Rq,, = 0.084 (0.1) 

which seems to justify our approach. 

R was introduced originally to take account of the anti-shielding caused by the overlap of charge distributions in 
the immediate neighbourhood of the target atom. 
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It should be noted that Royy = 0.081 and i?o,^ = 0.101 showing that R is not necessarily just a simple scalar 
parameter. However for the more symmetrical situated copper ions in the lattice with an axially symmetric EFG all 
components Rcua are identical. 

The apparent correspondence of our calculations with the terms of the Sternheimer equation to obtain R above is 
unsatisfactory in that the quantity 7 has had to be absorbed into the "target free ion", in other words, ^Vu = f{'y). 
Therefore i? is a function of 7, a result which questions the usefulness of either. This doubt is reinforced for 7 (see 
Ref.cj) since the large values necessary for the parameter clearly appear unsuitable to regard the anti-shielding as a 
perturbation. 

However before we disregard 7 entirely we have calculated it independently, assuming that the SE equation is 
correct, using first-principles calculations on a cluster used to model the effects of doping. The model and results are 
described in Appendix and agree with estimates in the literature that 7 is uncomfortably large. 

The evidence here points to the conclusion that the cuprate compounds cannot be satisfactorily analysed using 
formulae based on first-oder perturbation theory. 



C. EFG contributions from region I 

This section discusses the contributions to the EFG from the local, on-site term, ^M^^{Rj), first in a simplifying 
approach and then rigorously. This is important because V'°^°' is assumed to be "exact" in the Sternheimer equation. 

The basis functions at each centre are normally chosen to be radial functions, multiplied by a spherical harmonic, 
that is, hydrogen-like. The s-functions are spherical, the p- functions always occur as a group p^, Py and Pz, the 
d- functions as the group dz^-r^/^, d^x, dyz, d^^-y^, d^y, etc. 

A simplifying approach is based on the following argumentation. If all the functions in a group are equally occupied 
the associated electron density is spherical and the contribution to the EFG will be zero. So we are only interested in 
those functions which form part of the non-spherical density. So an electron configuration in the valence orbital 
will be analogous to the configuration or p^ (if both electrons have the same spin) that is, equivalent to a single 
hole in a spherical density. Of course here the ion will be in a crystal field and the degeneracies within each group 
may be lifted in which case we can be more definite than saying just p^ but, depending on the choice of axes, p^. 



This simphfies Eq. ( |3l8| ) to just ^Ml^{^{Rj) where 77 is the orbital which introduces the asymmetry into the spherical 



distribution about centre J. It is possible that the asymmetry is caused by the /3- spin orbital so the term we should 
use ^M^-'i^{Rj) as would be case in p^ for example, a (3 hole in the spherical distribution. 

Even ^M^^ ^{Rj) is simplified because in this discussion only one orbital on centre J is involved which we denote by 
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'Gt = 'Ml^{Rj) = 4%4%D'Ij^^{J). (5.4) 
This is therefore the required result subject to all these approximations introduced above 

V,,{Rj) = ^G^ (xng< >^ (5.5) 
where Ug is the occupancy of orbital g. In particular 



VxxiO) = -[2x n{2px) - n{2py) - n(2p,)] < r'^ (5.6) 



for the oxygen and 



2 

Vzz{Cu) ^ -[2x n(3rfj.2_j^2) + 2 X n{3dxy) 

- n{3dzx) - n{3dyz) - 2 X niid^z^^r-^)] < r^^ >3d (5.7) 



for the copper, and similarly for the other co mpon ents of Vi 



The right hand side of these equations ( p.q , |5.7| ) is just the form of the equation often used to estimate the occu- 
pancies, n, of orbitals. < >g is taken as the value for the atom or ion from calculation. This may be a poor 
estimate since < > is expected to be affected significantly by the crystal field which will tend to lift degeneracies 
and concentrate the asymmetries along particular directions. 

In the rigorous approach, the on-site matrix elements ^ are determined in the following way. The expectation 
values < >m,a are performed analytically. The coefficients c:^'^^ are given by the 7, a eigenvectors of the selfcon- 
sistent field equations. It must be emphasised that in the basis set 6-311G there are three radial functions for each of 
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the five 3d orbitals of Cu and three for each of the 2p orbitals. We remark that the relation 



< r 



>m,a— ^M^^^{Rj) 



which is valid for an orbital m showing 2p2;-symnietry implies the relation 



(5.8) 



(5.9) 



after performing the average over all the orbitals having the same symmetry. The resulting values for the oxygen in 
the above mentioned Cu9042 cluster are given in Tables III and E^. 



We first note that the expectation values < r~'^ > differ for the three p-orbitals and the spin projections by several 
percents in contrast to the assumptions in the simplified approach. It is evident that < >x differs from the other 
components since the bonding is along the x-direction, but also < r^'^ >y and < r^^ >z differ. This is due to a 
nonsymmetric distribution of the electron densities as it is plotted in Fig. |^ for the 2p electrons in oxygen. Furthermore 
the values < >~ 4 are about 10 % larger than those obtained from calculations on isolated atoms or ions. In 
Table IV we have also given the MuUiken partial charges (see Eq. ( 3.11| )) which show that the two are effectively 
fully occupied. It is only the 2p„ AO which is involved in the bonding and the convey of spin density from the copper 
to the ligand. 

The results of the analogous analysis of the contributions to the EFG and hyperfine dipole tensor for the Cu target 
ion are co llecte d in Tables ^ and The values called "remainder(s,d)" come from terms in the evaluation of matrix 
elements (3.17) where one basis function is s-like and the other d-like. The contributions from the Cu p-type orbitals 



to the EFG are substantial. This is due to the large < r~'^ > values. The occupancies of these orbitals, however, are 
close to one. The same applies to the three d orbitals with t2g symmetry as is seen from the MuUiken partial charges 
in Table VI. As expected, the distinguished AO is the 3d^2_y2 accompanied with some polarization in the Sd^^z^^,.^. 
Note, again, that < ~ 8, in contrast to the value ~ 6 used in approximate procedures. 

In a similar way the contributions ^^G^ and ^^Gg from region II can be analyzed. The results are given in 



Appendix ^. 



D. Electron- hole symmetry 



It should be emphasized that the occupations are determined by the expansion coefficients of the occupied 
MO into the individual AOs. To connect the results of the ab initio calculations with EFG analyses using the hole 
picture we can identify the unoccupied MOs which lie lowest in energy as contributions from "holes" . In particular, in 
our example of a CU9O42 cluster, there are nine such unoccupied MOs above the highest occupied molecular orbital 
(HOMO) which all show predominantly 3d^2_y2 and 2p^ character on the copper and oxygen rows, respectively. If 
we assume that they were occupied we wou ld ge t contributions t o the EFG and hyperfine tensor which we define by 



Vij and Ty. 



These are collected in Tables |VII| (for oxygen) and VIII (for copper) together with the values Vij and 



Tij as calculated from all occupied MOs (see Table 
For an isolated ion, one has the relations 



(5.10) 



and 



0. 



(5.11) 



For ions in t he cl uster, E q. ( 5.10) is n ot necessarily corre ct sin ce the environment is generall y non -spherical as is 
shown in Sec. VE. Tables VII and VIII show that relation ( |5.1l| ) approximately holds but that (5.10) is not fulfilled 
In this respect we remark that a calculation on the cluster (CuOg)^"'^^ (see Fig. |l|) where the Cu is nominally in a d 
state, yields Vz^ = -2.412. 
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E. Density matrix argument 



Since, contrary to the usual assumption, Eq. (5.10) is not correct under all conditions we will present the detailed 
theoretical background in this Section. We will not initially refer to the electron spin projection and so the terminology 
will be identical to that used in Sec. HID save for the a subscript. 
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So from Sec. |IIID| Eq. iK2^ 



(5.12) 



wher e c i s the matrix of MO coefficients collected in columns. The overlap matrix is defined as (using 0=1 in 
Eq. (p5|)) 



(5.13) 



However it is more convenient for us to use an orthogonal, but entirely equivalent, basis set of orbitals. We will label 
these by the column matrix B' such that 



B' = S-^B 



and 



=< B'IB't > 

where is the unit matrix of dimension Uc- The MO's are now written as 

$ = dB' 



(5.14) 



(5.15) 



(5.16) 



where d = cSa, the MO coefficients in the orthogonal basis. We now define spin MO's to hold electrons with 
spin projection and <i>^, to hold electrons with spin projection — i. Since each spin orbital can only hold one 
electron we need two density matrices, P and Q, to describe the a-spin and /3-spin densities respectively. 



P = d„I„dt 



Q 



d/jl/jdj. 



(5.17) 



where is a diagonal m x m matrix with I's for each occupied a-spin and zeroes otherwise. Ip is similar. The charge 
density matrix, necessary to calculate the EFGs, is given by P + Q and the spin density matrix, necessary to calculate 
the hyperfine tensor, is given by P — Q. In keeping with previous practice we can evaluate the "hole" density. The 
"hole" density is simply the total empty Hilbert space !,„ minus the P or Q. The charge "hole" density is 



and the spin "hole" density is 



Irn - Q = 2I„, - (P + Q) 



(Im-Q) = -(P-Q). 



(5.18) 



(5.19) 



Since the former includes the diagonal matrix 21^. this can make a contribution to the "hole" EFG calculation. If the 
Hilbert space is not spheri cal, o r at least does not possess cubic symmetry, the contribution will be non-zero. So the 
relation given in equation 5.10 is not strictly valid in ions where the degeneracy of the d-type orbitals is lifted. The 
degeneracy is only lifted if the symmetry is less than cubic. In all cases of practical interest 



(5.20) 



On the othe r hand the spin "hole" density is simply the negative of the spin density matrix leading to a verification 
of Eq. (ISTll) 



(5.21) 



We can use this argument to explain why, in Tables [VI| and VIII, the difference between the EFG's calculated from 
the occupied orbitals and the EFG's approximately calculated from selected unoccupied orbitals differ more markedly 
than the difference between equivalent calculations for the hyperfine tensors. If all the unoccupied orbitals are taken 
there is no difference between the occupied calculation and unoccupied calculation for the hyperfine tensor. 



VI. DISCUSSION AND CONCLUSIONS 



The number of problems which can be solved exactly by wave mechanics is very small and the perturbation method, 
originally devised for classical systems, was developed. Essentially, used in its less rigorous form, the problem is reduced 
to identifying that part which is well understood and treating the rest as a perturbation. In ionic crystals where the 
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properties, for example the EFG, of one of the ions (the target ion) is of interest, the purpose is to try and predict 
how changes to the properties of the isolated target ion can be accounted for by perturbation from its environment. 
In the mathematics of perturbation theory the changes to the target ion wave function can be achieved by mixing 
in the excited states of the target ion. Since the excited states form a complete set of functions this is always true 
although probably a very inefficient process. 

An environment of ions (point charges) contributes to the EFG at the target ion but will also interact with the 
electrons of the target ion to cause a distortion which in turn changes the electron contribution to the EFG. Such 
distortions should be easily simulated by judiciously mixing in the excited states of the target ion with its ground 
state. However if the possibility of covalent bonding occurs two problems seem to arise. Firstly the overlap of orbitals 
with nearest neighbour ions to the target ion and secondly the possibility, in a MuUiken population sense, of a transfer 
of electronic charge. Although in principle this can be accommodated by including the excited states of the target 
ion it is hardly a small perturbation questioning the applicability of the method to crystals where the possibility that 
covalent bonding occurs. 

The Sternheimer equation (Sec. II) uses a first-order perturbation theory type a rgume nt t o obt ain the individual 
terms which we attempt to correlate with different regions of the crystals (see Sees. [II B and HID ) . Star ting from a 
lattice of ionic charges the contribution to the EFG can be easily calculated at the target ion (Sec. IIIC). Of course 
the target ion, assumed to b e a point charge makes no contribution. The electronic structure about the target ion 
is very important (Sees. VD and VE) as long as it is not spherically symmetric. Since the crystal lattice interacts 
with the target ion (a crystal field) any asymmetry in the electron distribution (for example unfilled shells) will be 
significant since the crystal field will lift some degeneracies. The V^f""^ term is therefore crucial and fortunately is 
easily amenable to accurate calculation and transferable for the same ion to other crystals with a different chemical 
constitution. Although this crystal field could be guessed as being a small perturbation this is clearly not supported 
by the large value of 7 calculated here which are the same order of magnitude as those obtained "experimentally" 
(see also Appendix ^ . 

However the shielding parameter R (Sec. VB), intended to take account of "overlap" with nearest neighbour ions 
is rather more difficult to justify as a first order perturbation parameter. The "overlap" with the nearest neighbour 
orbitals could potentially lead to large electron density distortions, particularly of the outer shells, and also significant 
charge transfer. A better representation of this intuitive picture is the multi-center model which is at the core of 
molecular orbital theory. Unfortunately this complicates the interpretation of R and reduces its usefulness as does 
the conclusion that R and 7 are not independent. 

The use of experimental or theoretical data from isolated ions has long been a method of extracting information 
from a crystal system. Even without the complications of crystal fields or "overlap" the very existence of the lattice 
surrounding the target ion produces an unyielding restrictive cage from Pauli's exclusion principle. Ions can be 
attributed an ionic radius which apparently determines the structures of many ionic crystals. Any transfer of electronic 
charge onto the ion will hardly be able to use valence orbitals of the expected "free ion" size. This will no doubt 
contract the inner shell orbitals to compensate changing the experimental "free ion" EFGs. 

We conclude that the perturbed ion approach which results in the Sternheimer equation is inappropriate for cuprate 
crystals which are common in high-temperature superconducting materials mainly due to the significant covalent 
bonding in the Cu02 planes. However this has other consequences since the perturbation model also suggests an easy 
method to estimate the "holes" in the electron structure from the EFGs, whose distribution in turn is essential for 
models of superconductivity itself. We have shown that these estimates are probably wrong and at the very least 
their values should be reassessed. Therefore precise information on the charge and spin density distributions in copper 
oxides is necessary and EFGs, determined by nuclear quadrupole resonance spectroscopy, can help to provide this 
information. It is necessary, however, that they are analysed in a more sophisticated manner than with point-charge 
models. 
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APPENDIX A: CONTACT DENSITIES 



The regional partitioning for the evahiation of the EFG tensors and the dipolar hyperfine tensors applies also for 
the contact interaction. The corresponding results are given in this appendix. 
We denote the contact density for the target nucleus J as 




where the sum extends over the occupied MOs and perform the same regional partitioning as in Sec. IV. The total 



contributions to I?„s(Cu) and D„s(0) for the different s-like AOs are listed in Table [X with the small contributions 



from regions II and III given in parentheses. Since the expectation values | '4'ns{Rj) I have nearly the same values 
for spin up and down projections we can describe the results also in terms of partial polarizations fns according to 

Dns{Rj) - Y I AM) ? fns- (A2) 

Note that these results refer to maximal spin-multiplicity. Thus, the values for D(C\x) and D{0) include the 
transferred hyperfine fields from the four and two ncc^p^^ neighbour copper ions, respectively. These transferred 
hyperfine fields have been discussed extensively in Refs 



APPENDIX B: ESTIMATE OF 7 PARAMETER 

We have performed several cluster calculations where point charges q have been added to the La'^+ pseudopoten- 
tials at positions X, Y, and Z for the clusters Cu06/Cu4Laio and Cu20ii/Cu6Lai6 (see Fig. ||). For the cluster 
Cu06/Cu4Laio the target ions are the central Cu and the planar O on the x-axis whereas for Cu20ii/Cu6Lai6 the 
target ions are the Cu to the right and the central O. 

Since these additional charges are in region III, the differences in the calculated EFG tensors are then identified 
with the term /Wu^q) = (1 - 7)AM/(qt) 



The results are collected in Tables PCII| and XIII , It is seen that these 7 values are unreasonably large and that these 
"lattice" contributions in the SE cannot be used at all. What really happens is that the additional charges distort 
and polarize the nearby ions (oxygens in the present cases) which in turn then influence the target ion. 



APPENDIX C: CONTRIBUTIONS FROM REGION II 



For completeness we collect here the results of the analysis of the contributions from region II (see Sec. V C ). With 
the oxygen as target nucleus, the values of ^^G^^ and ^^G^' are given in Table Note that the contributions assigned 
to s-character are due to matrix elements of the operator V between s-type functions centred at the oxygen and d-type 
functions centred at the neighboring copper nuclei. 



In Table XI the contributions from region II for the copper target nucleus are given. 
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TABLE I: Diagonal elements of the tensors G from the contributions of on-site (I), on-site/off-site (II) and off-site AOs (III) 
for spin projections a and (3 for Cu and O. Contributions from nuclear charges W , EFG tensor V and hyperfine tensor T . 







Cu 











G 


XX 


yy 


zz 


XX 


yy 


zz 




0.432 


0.432 


-0.864 


-0.141 


0.130 


0.011 




-1.266 


-1.266 


2.532 


-0.787 


0.436 


0.351 




0.073 


0.073 


-0.146 


0.018 


-0.013 


-0.005 




0.102 


0.102 


-0.204 


0.060 


-0.033 


-0.027 


111,^ 


0.278 


0.278 


-0.556 


1.367 


-0.418 


-0.949 


III^ 

Gg 


0.262 


0.262 


-0.524 


1.286 


-0.390 


-0.896 


w 


-0.522 


-0.522 


i.Oii 


-2. ()!):! 


l).8:i2 


1.861 


V 


-0.642 


-0.642 


1.283 


-0.890 


0.545 


0.345 


T 


1.685 


1.685 


-3.370 


0.685 


-0.314 


-0.371 



TABLE II: Contributions to the EFG and the hyperfine coupling tensor from the different regional partitions. 



Region V"..(Cu) r..(Cu) 14.x (O) Vyv{0) V;.(Q) T^^Q) r^^(Q) r,.(Q) 
I 1.668 -3.396 -0.928 0.566 0.362 0.646 -0.306 -0.340 
II -0.350 0.058 0.078 -0.046 -0.032 -0.042 0.020 0.022 
III + W -0.036 -0.032 -0.040 0.024 0.016 0.081 -0.028 -0.053 
Total 1.282 -3.370 -0.890 0.544 0.346 0.685 -0.314 -0.371 



TABLE III: Contributions of on-site AOs (region I) for spin projection a and /J for the planar oxygen. 







a 












XX 


yy 


zz 


XX 


yy 


zz 


Px 


3.046 


-1.523 


-1.523 


2.370 


-1.185 


-1.185 


Py 


-1.613 


3.226 


-1.613 


-1.592 


3.184 


-1.592 


Pz 


-1.574 


-1.573 


3.147 


-1.564 


-1.564 


3.128 


'G 


-0.141 


0.130 


0.011 


-0.787 


0.436 


0.351 
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TABLE IV: Occupations ^N{2pj) and MuUiken partial charges p of the 2p orbitals, and averaged values of < r ^ >. 







2px 


2py 


2pz 






0.921 


1.018 


1.006 






0.737 


1.013 


1.006 


Pr. 




0.914 


0.999 


1.001 


Pd 




().7cs:i 


0.991 


i.OOl 






4.135 


3.961 


3.911 




>P 


4.044 


3.933 


3.885 



TABLE V: Contributions of on-site AOs (region I) for spin projections a and (3 for the central copper. 





a 


P 


Remainder(d,s) 


0.214 


0.232 


P 


-0.651 


-0.599 




-4.839 


-1.325 


d^2_^2/3 


4.453 


4.277 


dxy 


-4.563 


-4.538 


dzx^ dyz 


4.522 


4.485 


^Gzz 


-0.864 


2.532 



TABLE VII: Electric field gradients Vu and the hyperfine tensors Tu for the occupied states and Vu and Tu from the unoccupied 
states close but above the Fermi energy for the planar oxygen atom. 



i 


Vu 


F« 




T a 


X 


-0.873 


0.624 


0.685 


-0.624 


y 


0.563 


-0.284 


-0.314 


0.284 


z 


0.310 


-0.340 


-0.371 


0.340 



TABLE VIII: Electric field gradients Vu and the hyperfine tensors Tu for the occupied states and Vu and Tu from the 
unoccupied states close but above the Fermi energy for the central copper nucleus. 



i 


Vu 


Vu 




T a 


x,V 


-0.623 


1.670 


1.685 


-1.670 


z 


1.246 


-3.340 


-3.370 


3.340 



TABLE VI: Occupations and MuUiken partial charges p of the 3d orbitals, and averaged values of < r ^ >. 







3dj.2_y2 


3(^2,2 _^2/3 






Sdyz 






1.030 


0.966 


1.000 


0.998 


0.998 






0.297 


0.939 


1.000 


0.998 


0.998 


pa 




0.999 


0.971 


0.997 0.997 


0.997 


P0 




0.358 


0.946 


0.997 


0.997 


0.997 




>a 


8.224 


8.066 


7.982 


7.926 


7.926 


■<r-3 


>/3 


7.821 


7.976 


7.939 


7.864 7.864 
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TABLE IX: Expectation values of the s-like AOs, |7/;„3|^, contact densities D„s, and polarizations f„s at Cu and O, respectively. 
For D„s the total values are given with the contributions from regions II and III in parentheses. 



n 


V'ns(Cu) 


1" Dn.(Cu) 


fns[%] |V-(0)|" 


Dns(O) 


fns[%] 


1 


7300 


-0.053 (0.000) 


-8.66 X 10"^ 141 


-0.262 (0.001) 


-0.0222 


2 


725 


-3.637 (0.002) 


-5.99 X 10"^ 6.77 


1.516 (0.055) 


2.673 


3 


107 


2.575 (0.013) 


0.288 






4 


2.35 


2.037 (-0.017) 


10.35 







TABLE X: Contributions to "G from region II for spin projection a and /? for the planar oxygen. 







a 












XX 


yy 


zz 


XX 


yy 


zz 


s 


-0.004 


0.001 


0.003 


0.002 


-0.002 


-0.001 


Px 


0.021 


-0.011 


-0.010 


0.057 


-0.028 


-0.028 


Py 


0.001 


-0.004 


0.003 


0.001 


-0.004 


0.003 


Pz 


0.000 


0.001 


-0.001 


0.000 


0.001 


-0.001 




0.018 


-0.013 


-0.005 


0.060 


-0.033 


-0.027 



TABLE XI: Contributions to "G from region II for spin projections a and /? for the central copper. 





a 


P 


Remainder 


0.000 


0.000 


s 


0.001 


0.003 


P 


0.011 


0.011 




0.013 


-0.025 


C^22-r2/3 


-0.172 


-0.194 




0.001 


0.001 


dzxi dyz 


0.000 


0.000 


*^zz 


-0.146 


-0.204 



TABLE XII: Calculated values of 7 for the cluster in Fig. h|(a) with additional point charges at positions X and Y. 



X 


Y 7(Gu) 7(0) 


-0.1 


-0.1 -51.1 +19 


-0.1 


-51.3 +20 


-1.0 


-50.9 +19 


-1.0 


-1.0 -52.1 +23 



TABLE XIII: Calculated values of 7 for the cluster in Fig. m(b) with additional point charges at positions X, Y, and Z. 



X 


Y 


Z 7(Git) 


7(0) 


-0.1 


-0.1 


-0.1 -37 


-33 





-0.1 


-0.1 -40 


-32 








-1.0 -26 


-33 



18 



(a) 



(b) 



•X 




FIG. 1: The Cu06/Cu4Laio and Cu2 0ii/Cu6Lai6 clusters. The notations X, Y, and Z are for later reference (see Sec. 
Appendix [B|) . 




FIG. 2: Illustration of the contributions I to III to expectation values in a Cu04 cluster. The Jth atom is the central Cu, 
whereas the neighboring 0-atoms denote the atoms K and L. The full curve limits the d-electrons of the central Cu and the 
dotted curves enclose the oxygen p-electrons. 



% O'-' ^ -V - 



- o 



Cu 



FIG. 3: Highest occupied molecular orbital for the Cu9O42/Cui2La50 cluster. 
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FIG. 4: Density distribution of the 2p oxygen electrons in a yi-plane perpendicular to the Cu-O-Cu connection line and 
through the oxygen atom. The arrows point along increasing densities. The equidensity lines close the density max;ima show 
that these maxima are larger in the y than in the z direction. 



